﻿using LxBasic;
using System.Diagnostics;

namespace LxCore
{
    public class SIMPLEPressure : MatrixBase
    {
        //local
        public int ReferenceCell;//如果没有压力边界，会导致缺少一个约束，则需要设置一个参考压力点
        public double ReferencePressure;//记录下来，避免多次改设置偏了
        public DiagonalMatrixR3C3[] InnDc;//Dc，见591。用以求解Df,DDf。DDf见P596//Dc在保存的时候，剔除松弛的影响！
        public Vector3d[] Grad_LastStep;
        public PhysicalFace[] InnPhyFaces;
        public PhysicalFace[] BouPhyFaces;
        //share
        public SIMPLEFields InSimpleFields { get; set; }
        public double[] InnX_LastStep;
        public double[] BouX_LastStep;
        public override void CheckCreate(FvZone zone)
        {
            base.CheckCreate(zone);
            CheckCreator.ArrayStrict(ref InnDc, InnCells.Length);
            CheckCreator.ArrayStrict(ref InnPhyFaces, InnFaces.Length);
            CheckCreator.ArrayStrict(ref BouPhyFaces, BouFaces.Length);
            CheckCreator.ArrayStrict(ref Grad_LastStep, InnCells.Length);
        }
        public override bool IsSymmetric()
        {
            return false;//设置参考压力将导致非对称矩阵
        }

        public virtual void SetZone(IField<double> iniValue, IParameter<double> relaxFactor)
        {
            this.IniValue = iniValue;
            this.RelaxFactor = relaxFactor;
        }
        public override void Initialize()
        {
            if (!IsXInit)
            {
                Vector3d center;
                MultiThreads.RunParallelTask(ini1_IniX_CalCenter);
                {
                    (double sum1, double sum2, double sum3) = ThreadsMath.MemoryRecallClear();
                    center = new Vector3d(sum1, sum2, sum3) * (1 / Zone.Inner.Volumn);
                }
                double MinLength = double.MaxValue;
                MultiThreads.RunParallelTask(ini2_Grad_FindRefCell);
                IsXInit = true;
                IsGradInit = true;
                void ini1_IniX_CalCenter(ThreadInfo threadInfo)
                {
                    SetXs(IniValue, Zone.Inner, threadInfo);
                    var mbs = InSimpleFields.SimpleBounds;
                    for (int i1 = 0; i1 < mbs.Length; i1++)
                    {
                        switch (mbs[i1])
                        {
                            case SIMPLEFields.PressureInOut pb:
                                SetXs(pb.Pressure, Zone.Bounds[i1], threadInfo);
                                break;
                            default:
                                SetXs(IniValue, Zone.Bounds[i1], threadInfo);
                                break;
                        }
                    }
                    //为了求中心，先求体积*坐标之和
                    {
                        (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
                        Float128 s1 = new Float128(), s2 = new Float128(), s3 = new Float128();
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnCell cell = ref InnCells[i1];
                            Vector3d v = cell.Volumn * cell.Center;
                            s1 += (v.X); s2 += (v.Y); s3 += (v.Z);
                        }
                        ThreadsMath.MemoryAdd(s1, s2, s3);
                    }
                }
                void ini2_Grad_FindRefCell(ThreadInfo threadInfo)//和上面计算面积的要对应
                {
                    CalGradient_LSQ(threadInfo);
                    {
                        double minLength = double.MaxValue;
                        int minLengthCell = 0;
                        (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
                        for (; i1 < endp1; i1++)
                        {
                            double length = (InnCells[i1].Center - center).LengthSquared();//实际比较的长度平方
                            if (length < minLength)
                            {
                                minLength = length;
                                minLengthCell = i1;
                            }
                        }
                        lock (MultiThreads.GlobalLock)
                        {
                            if (minLength < MinLength)
                            {
                                MinLength = minLength;
                                ReferenceCell = minLengthCell;
                                ReferencePressure = InnX[ReferenceCell];
                            }
                        }
                    }
                }
            }
            if (!IsGradInit)
            {
                MultiThreads.RunParallelTask(CalGradient_LSQ);
                IsGradInit = true;
            }
        }
        public void SetReferencePoint()
        {
            int cellIndex = ReferenceCell;
            ref FvInnCell cell = ref InnCells[cellIndex];
            int i1 = 0;
            for (; i1 < cell.InnOFaceEnd; i1++)
            {
                AInnO[cell.AdjFaces[i1]] = 0;
            }
            for (; i1 < cell.InnNFaceEnd; i1++)
            {
                AInnN[cell.AdjFaces[i1]] = 0;
            }
            for (; i1 < cell.AdjFaces.Length; i1++)
            {
                ABouH[cell.AdjFaces[i1]] = 0;//如果cell位于交界面附近，BouAH可能会被赋值，出现错误
            }
            B[cellIndex] = ReferencePressure * ADiag[cellIndex];
        }
        public override void FormABOnFaces(ThreadInfo threadInfo)
        {
            {
                (int i1, int endp1) = Zone.Inner.FaceIdx.Split(threadInfo);
                var density = InSimpleFields.Density.GetFaceVarArr(Zone.Inner);
                for (; i1 < endp1; i1++)
                {
                    ref FvInnFace face = ref InnFaces[i1];
                    ref PhysicalFace mface = ref InnPhyFaces[i1];
                    double den = density.GetFaceVar(i1);
                    var Df = den * face.CalAverage(InnDc);
                    mface.SetFace(Df, face);
                    AInnO[i1] = AInnN[i1] = -mface.gDiff;
                    InnFaceShare[i1] = mface.gDiff * (InnX_LastStep[face.OSideCell] - InnX_LastStep[face.NSideCell]);
                }
            }
            var mbs = InSimpleFields.SimpleBounds;
            for (int i1 = 0; i1 < mbs.Length; i1++)
            {
                switch (mbs[i1])
                {
                    case SIMPLEFields.PressureInOut pb:
                        set_ALU_PressureBound(Zone.Bounds[i1], pb.Pressure, threadInfo);
                        break;
                    case SIMPLEFields.NoSlipWall noslip:
                        setFactorOfMassBound(Zone.Bounds[i1], noslip.FactorsForPressure, threadInfo);
                        break;
                    case SIMPLEFields.Symmetry sym:
                        setFactorOfMassBound(Zone.Bounds[i1], sym.FactorsForPressure, threadInfo);
                        break;
                    case SIMPLEFields.VelocityInlet velInlet:
                        setFactorOfMassBound(Zone.Bounds[i1], velInlet.FactorsForPressure, threadInfo);
                        break;
                }
            }
            void set_ALU_PressureBound(FvBound bound, IField<double> pressure, ThreadInfo threadInfo)
            {
                (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
                var density = InSimpleFields.Density.GetFaceVarArr(bound);
                var p = pressure.GetCellVarArr(bound);
                double[] BouMassFlux = InSimpleFields.MassFlux.BouMassFlux;
                var bouVel = InSimpleFields.Velocity.GetBouCellVelVars();
                for (; i1 < endp1; i1++)
                {
                    ref var face = ref BouFaces[i1];
                    ref PhysicalFace mface = ref BouPhyFaces[i1];
                    double den = density.GetFaceVar(i1);
                    var Df_average = den * InnDc[face.OSideCell];
                    mface.SetFace(Df_average, face);
                    double pbou = 0;
                    double massflux = BouMassFlux[i1];
                    if (double.IsNegative(massflux))
                    {//流入
                        pbou = p.GetCellVar(i1) - 0.5 * den * bouVel.GetCellVar(i1).LengthSquared();
                    }
                    else
                    {//流出
                        pbou = p.GetCellVar(i1);
                    }
                    BouX[i1] = pbou;// = p.GetCellVar(i1);
                    ABouH[i1] = -mface.gDiff;
                    BouFaceShare[i1] = mface.gDiff * (pbou + InnX_LastStep[face.OSideCell] - BouX_LastStep[i1]);
                }
            }
            void setFactorOfMassBound(FvBound bound, Vector2d[] factors, ThreadInfo threadInfo)
            {
                FormALUb_0(bound, threadInfo);
                (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);             
                double[] bouP = BouX; double[] innP = InnX; Vector3d[] Grad = this.Grad;
                double[] bouFlux = InSimpleFields.MassFlux.BouMassFlux;
                var innVel = InSimpleFields.Velocity.GetInnCellVelVars();
                var density = InSimpleFields.Density.GetFaceVarArr(bound);
                int start = bound.FaceIdx.Start;
                for (; i1 < endp1; i1++)
                {
                    ref var face = ref BouFaces[i1];
                    ref var pFace = ref BouPhyFaces[i1];
                    double den = density.GetFaceVar(i1);
                    var Dc = InnDc[face.OSideCell];
                    pFace.SetFace(den * Dc, face);
                    double flux_b = bouFlux[i1] - den * Vector3d.DotProduct(innVel.GetCellVar(face.OSideCell) + Matrix3dMath.Product(Dc, Grad_LastStep[face.OSideCell]), face.Normal);
                    factors[i1 - start] = new Vector2d(
                        -pFace.CalNonOrthogonalityCorrectionOfBound(ref face, Grad, innP, bouP[i1]) / pFace.gDiff,
                        -flux_b / pFace.gDiff
                        );
                    //BouFaceShare[i1] = flux_b;//与热流密度边界不一样，不存在那样的关系，没研究透彻。。。
                }
            }
        }
        public void UpdateABOnFaces(ThreadInfo threadInfo)
        {
            {
                (int i1, int endp1) = Zone.Inner.FaceIdx.Split(threadInfo);
                for (; i1 < endp1; i1++)
                {
                    ref FvInnFace face = ref InnFaces[i1];
                    ref PhysicalFace mface = ref InnPhyFaces[i1];
                    InnFaceShare[i1] = mface.CalNonOrthogonalityCorrectionOfInner(ref face, Grad, InnX)
                        - Vector3d.DotProduct(mface.Normal, face.CalCFCorrectedGradientOfInner(Grad_LastStep, InnX_LastStep));
                }
            }
            var mbs = InSimpleFields.SimpleBounds;
            for (int i1 = 0; i1 < mbs.Length; i1++)
            {
                switch (mbs[i1])
                {
                    case SIMPLEFields.PressureInOut pb:
                        set_ALU_PressureBound(Zone.Bounds[i1], pb.Pressure, threadInfo);
                        break;
                    case SIMPLEFields.NoSlipWall noslip:
                        setFactorOfMassBound(Zone.Bounds[i1], noslip.FactorsForPressure, threadInfo);
                        break;
                    case SIMPLEFields.Symmetry sym:
                        setFactorOfMassBound(Zone.Bounds[i1], sym.FactorsForPressure, threadInfo);
                        break;
                    case SIMPLEFields.VelocityInlet velInlet:
                        setFactorOfMassBound(Zone.Bounds[i1], velInlet.FactorsForPressure, threadInfo);
                        break;
                }
            }
            void set_ALU_PressureBound(FvBound bound, IField<double> pressure, ThreadInfo threadInfo)
            {
                (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
                for (; i1 < endp1; i1++)
                {
                    ref var face = ref BouFaces[i1];
                    ref PhysicalFace mface = ref BouPhyFaces[i1];
                    double pbou = BouX[i1];
                    BouFaceShare[i1] = mface.gDiff * pbou + mface.CalNonOrthogonalityCorrectionOfBound(ref face, Grad, InnX, pbou)
                        - Vector3d.DotProduct(mface.Normal, face.CalCFCorrectedGradientOfBound(Grad_LastStep, InnX_LastStep, BouX_LastStep[i1]));
                }
            }

            void setFactorOfMassBound(FvBound bound, Vector2d[] factors, ThreadInfo threadInfo)
            {
                (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
                FormALUb_0(bound, threadInfo);
                double[] bouP = BouX; double[] innP = InnX; Vector3d[] gradP = this.Grad;
                int start = bound.FaceIdx.Start;
                for (; i1 < endp1; i1++)
                {
                    ref var face = ref BouFaces[i1];
                    ref var pFace = ref BouPhyFaces[i1];
                    factors[i1 - start].X = -pFace.CalNonOrthogonalityCorrectionOfBound(ref face, gradP, innP, bouP[i1]) / pFace.gDiff;
                    //BouFaceShare[i1] = -factors[i1 - start].Y * pFace.gDiff;
                }
            }
        }

        public override void FormABOnCells(ThreadInfo threadInfo)
        {
            double[] massImb = InSimpleFields.MassFlux.MassImbalance;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                ref FvInnCell cell = ref InnCells[i1];
                ADiag[i1] = CalAdiag(ref cell);
                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) - massImb[i1];
            }
        }
        public void UpdateABOnCells(ThreadInfo threadInfo)
        {
            double[] massImb = InSimpleFields.MassFlux.MassImbalance;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                ref FvInnCell cell = ref InnCells[i1];
                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) - massImb[i1];
            }
        }
        public void SaveDc(ThreadInfo threadInfo)
        {
            var Dc = InnDc;
            var U_ADiag = InSimpleFields.Velocity.U.ADiag;
            var V_ADiag = InSimpleFields.Velocity.V.ADiag;
            var W_ADiag = InSimpleFields.Velocity.W.ADiag;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                ref var cell = ref InnCells[i1];
                double volumn = cell.Volumn;
                Dc[i1] = new DiagonalMatrixR3C3(volumn / U_ADiag[i1], volumn / V_ADiag[i1], volumn / W_ADiag[i1]);
            }
        }
        public void SaveXAndGradLastStep(ThreadInfo threadInfo)
        {
            LxTool.CopyArray(InnX, InnX_LastStep, threadInfo);
            LxTool.CopyArray(BouX, BouX_LastStep, threadInfo);
            LxTool.CopyArray(Grad, Grad_LastStep, threadInfo);
        }
        public override void UpdateBoundsValue(ThreadInfo threadInfo)
        {
            var mbs = InSimpleFields.SimpleBounds;
            for (int i1 = 0; i1 < mbs.Length; i1++)
            {
                switch (mbs[i1])
                {
                    case SIMPLEFields.PressureInOut:
                        break;
                    case SIMPLEFields.NoSlipWall noslip:
                        predictVelBoundPressure(Zone.Bounds[i1], noslip.FactorsForPressure, threadInfo); break;
                    case SIMPLEFields.Symmetry sym:
                        predictVelBoundPressure(Zone.Bounds[i1], sym.FactorsForPressure, threadInfo); break;
                    case SIMPLEFields.VelocityInlet velInlet:
                        predictVelBoundPressure(Zone.Bounds[i1], velInlet.FactorsForPressure, threadInfo); break;
                }
            }
            void predictVelBoundPressure(FvBound bound, Vector2d[] factors, ThreadInfo threadInfo)
            {
                int start = bound.FaceIdx.Start;
                (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
                for (; i1 < endp1; i1++)
                {
                    ref var face = ref BouFaces[i1];
                    ref Vector2d factor = ref factors[i1 - start];
                    BouX[i1] = InnX[face.OSideCell] + factor.X + factor.Y;
                }
            }
        }
        public void CalGradAndCorrect_SIMPLE_1(ThreadInfo threadInfo)
        {
            var velo = InSimpleFields.Velocity;
            var U_x = velo.U.InnX;
            var V_x = velo.V.InnX;
            var W_x = velo.W.InnX;
            var relax = RelaxFactor.GetPara();
            var relax_1 = 1 - relax;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                var grad = InnCells[i1].CalGradient_LeastSquare(InnFaces, BouFaces, InnX, BouX);
                Vector3d gradp = grad - Grad_LastStep[i1];
                Vector3d _velp = Matrix3dMath.Product(InnDc[i1], gradp);
                U_x[i1] -= _velp.X; V_x[i1] -= _velp.Y; W_x[i1] -= _velp.Z;
                Grad[i1] = grad - relax_1 * gradp;
                InnFaceShare[i1] = relax * InnX[i1] + relax_1 * InnX_LastStep[i1];
            }
            (i1, endp1) = new IndexRange(BouFaces.Length).Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                BouFaceShare[i1] = relax * BouX[i1] + relax_1 * BouX_LastStep[i1];
            }
        }
        public void CalGradAndCorrect_SIMPLE_2(ThreadInfo threadInfo)
        {
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            Array.Copy(InnFaceShare, i1, InnX, i1, endp1 - i1);
            (i1, endp1) = new IndexRange(BouFaces.Length).Split(threadInfo);
            Array.Copy(BouFaceShare, i1, BouX, i1, endp1 - i1);
        }
        public void CalGradAndCorrect_PISO(ThreadInfo threadInfo)
        {
            var velo = InSimpleFields.Velocity;
            var U_x = velo.U.InnX;
            var V_x = velo.V.InnX;
            var W_x = velo.W.InnX;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                var grad = InnCells[i1].CalGradient_LeastSquare(InnFaces, BouFaces, InnX, BouX);
                Vector3d gradp = grad - Grad_LastStep[i1];
                Vector3d _velp = Matrix3dMath.Product(InnDc[i1], gradp);
                U_x[i1] -= _velp.X; V_x[i1] -= _velp.Y; W_x[i1] -= _velp.Z;
                Grad[i1] = grad;
            }
        }
        public override void FromJsonKV(string key, IJsonReader jsReader)
        {
            switch (key)
            {
                case nameof(ReferenceCell): ReferenceCell = jsReader.ReadInt(); break;
                case nameof(ReferencePressure): ReferencePressure = jsReader.ReadDouble(); break;
                default: base.FromJsonKV(key, jsReader); break;
            }
        }
        public override void ToJsonKV(IJsonWriter jsWriter)
        {
            base.ToJsonKV(jsWriter);
            jsWriter.WriteInt(nameof(ReferenceCell), ReferenceCell);
            jsWriter.WriteDouble(nameof(ReferencePressure), ReferencePressure);
        }
    }
}
